Canonical and Microcanonical Ensemble Approaches to Bose-Einstein Condensation: 
The Thermodynamics of Particles in Harmonic Traps 
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The thermodynamic properties of bosons moving in a harmonic trap in an arbitrary number of 
dimensions are investigated in the grand canonical, canonical and microcanonical ensembles by ap- 
plying combinatorial techniques developed earlier in statistical nuclear fragmentation models. Ther- 
modynamic functions such as the energy and specific heat are computed exactly in these ensembles. 
The occupation of the ground or condensed state is also obtained exactly, and signals clearly the 
phase transition. The application of these techniques to fermionic systems is also briefly discussed. 



I. INTRODUCTION 

Recently there has been a renewed interest in the ther- 
modynamics of mesoscopic systems. In particular, the re- 
cent experimental observation of condensation phenom- 
ena in nanokelvin aggregates of atoms moving in mag- 
netic traps Q-Q] has encouraged a number of recent the- 
oretical investigations of this phenomenon extend- 
ing earlier investigations |^,|^]. These works assume that 
the systems can be adequately described by the grand 
canonical ensemble in which both the energy and particle 
number are allowed to fluctuate. The canonical ensemble 
is probably more representative of the experimental con- 
ditions and has also been recently investigated ]To| , |TT] | . 

In this work, we analyze the thermodynamics of bosons 
in the the grand canonical, canonical and microcanonical 
ensembles by applying particularly effective techniques 
that two of us have used before in describing nuclear 
systems statistically |l2|-|l4|]. This reproduces the re- 
cursive formula for the canonical partition function from 
a combinatorial argument which proffers several advan- 
tages over the earlier purely functional treatment [ ^0| . 
In particular, this interpretation enables a simple and 
exact determination of the expected occupancies of the 
energy levels in any of the ensembles. This allows the 
usual thermodynamic functions to be expressed exactly 
and succinctly, not just the partition function. Addition- 
ally, the method is related to the permutational decom- 
position of the liquid 4 He partition function proposed by 
Feynman |l^,|l6|. Previously, we have used this corre- 
spondence to specify partition functions for isotopic nu- 
clear systems |T7fl . In this paper, the correspondence 
enables as a new result an exact determination of the 
microcanonical partition function. 



II. THE PARTITION FUNCTION FOR 
IDENTICAL BOSONS OR FERMIONS 

In statistical mechanics the partition function entirely 
specifies the thermodynamics, so our first priority is to 
determine a method of computing the partition function 
in the various ensembles. We begin by considering the 



grand canonical partition function, and develop a com- 
binatorial method of extracting the canonical and micro- 
canonical partition functions from the larger ensemble. 
In the process we develop exact expressions for the ex- 
pected occupation of the energy levels in all three ensem- 
bles. 



A. Grand canonical partition function 

In the grand canonical ensemble, where the particle 
number and energy are unconstrained, the partition func- 
tion Z for a system of bosons occupying energy levels Ek 
with degeneracy gu is given by 
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where /x is the chemical potential. Let us think of this as 
a sum either over canonical Z n or microcanonical Z n ^E 
partition functions, i.e. 



Z = ^z n Z n = ^z n w E Z, l 



n,E , 
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where the fugacity z is given by z = exp{/i/fcsT} and 
w = cxp{ — 1/ksT}. Expanding Eq. (Q) as a power series 
in z,w yields 
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where we have introduced Xjk = 9j/k and Xk = 
Y]j w k£j Xjk for reasons which will become clear shortly. 
If we exponentiate and Taylor expand Eq. (||) we see that 

z - 1^ —\ <L : x 



7T01 



7T01=0 " 7I"02=0 

(z k W kej XjkY ]k 



7T02! 



En 



TTik 



(4) 



1 



where irjk > 0. 

Notice that the exponents of w, z in Eq. (^) are simple 
functions of {iTjk}- Since the exponent of w in the grand 
canonical partition function is the total energy E in the 
particular microcanonical component, the exponent of w, 
Sfc kEjTTjk should be equal to the energy E for that par- 
ticular state {ftjk}- This E — Ylj n 3 £ ji where rij is the 
number of particles in energy state j, so the state 7? is 
related to ft by 
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Combinatorially speaking, njk counts collections of k par- 
ticles in the jth energy state. In fact, {itjk '■ j fixed} is 
simply a cycle class decomposition of a permutation of 
the rij particles in the jth energy level. To see this, con- 
sider the part of the grand canonical partition function 
due to energy level j, namely the number of states of a 
system with rij bosons in a level with degeneracy gj, 



rij + gj - 1 
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Here we see that ( nj+ ^ j ) can be described in two ways. 
The first is just the isolated contribution from the above 
grand canonical partition function. The second is a sum 
over permutations of the rij particles, where ftjk(p) is 
the number of fc-cycles in the permutation p of rij . The 
equivalence of these two expression is due to Cauchy and 
Sylvester ||. 

Feynman considered such permutations in the conden- 
sation of liquid 4 Hc [[[TJ] and was able to obtain a form 
for the partition function by expanding the symmetrized 
density matrix as a sum over such permutations. Specif- 
ically, the sum iik = Y]j iTjk, the total number of cycles 
of length fc, figures prominently in Feynman's approach. 
Such permutational decomposing of the partition func- 
tion has also been useful in our research on isotopic nu- 
clear systems |Tij| , where the dual requirement of proton 
and neutron conservation was considered as restricting 
the permutational states to have been colored in a par- 
ticular way. Here we can accomplish a similar interpre- 
tation by coloring the clusters which form according to 
the energy levels they reside in. Parallels of permutation 
problems with cluster yields were also noted in [p^|Jl4[ . 

Since rtjk has a combinatorial interpretation and is re- 
lated to rij, it is important to obtain expectation values 
of it. This is possible since each term of Eq. (|j) can be 
interpreted as an unnormalized probability for the state 
{TTjfc}, the normalization being simply the partition func- 
tion. Equation (Q) then suggests we can evaluate such 
expectation values of irjk by taking various derivatives of 
the partition function. More specifically, it is true that 
(ir jk ) = (xjk/Z)(dZ/dx jk ) and (7r ife (7r /m - S jk ,i m )) = 



(xjkXim/Z)(d 2 Z/dxjkdxi m ), where (•} denotes a grand 
canonical ensemble average. Now according to Eq. (0), 
Z = exp V - fe z k w h£j Xjk, so that 
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We then see that 



(wjk) = z k W ks * Xj k , (8) 



(njk(7Tlm ~ Sjkjm)) = Z k+m W k ^ + m ^X jk 



(9) 



We can now "solve" the grand canonical partition func- 
tion now by using the combinatorial interpretation to 
specify w, z, i.e. by applying the constraints n — 
E,-K) = E,-fc kfak), E = Y,Anj)ej = Y, jk k(nj k )ej. 



The analysis illustrated above for a system of bosons 
applies equally well to a system of fermions. In this case, 
the appropriate partition function is 



InZ = Y^9j hi(l + eb-*iV kaT ) . 



(10) 



Expanding this yields Eq. (g) with Xjk = (—l) + gj/k. 
The fact that Xjk can be negative means that the combi- 
natorial interpretation is somewhat suspect, since some 
contributions to the partition function are negative. Nev- 
ertheless, the mathematics is sound and expectation val- 
ues of {rij} can be computed correctly by the same tech- 
nique used above. This fermionic case will be discussed 
in detail in another paper. 



B. Canonical partition function 

We can obtain the canonical partition function from 
the grand canonical partition function by making restric- 
tions on the sum appearing in Eq. (0), i.e. if we restrict 



the sum to terms where 



- zZji 



knjk we will obtain the 



canonical partition function. We could simplify some of 
the computations by defining irk = Yjftjk, so that the 
states are described by {rrk} where J2 k kirk — n. Mathe- 
matically speaking, such vectors {irk} describe partitions 
of the integer n |18| , where ~Kk is the number of time k 
appears in such a partition. However, sometimes it is 
more convenient to work instead with iXjk- 

Let us define tt(ji) = {iXjk '■ Sjfe ^ 7r j'fe = n }- Then the 
canonical partition function is given by 
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Since JZj^ kiTjk — n for each state in the ensemble, it is 



true that ^ - fe k{-Kjk)n = n, where (•)„ denotes a canon- 



ical ensemble average. The canonical ensemble average 
for TTjk can be computed by taking a derivative of the 
partition function (rrjk)n = (%jk/Z n )(dZ n /dxjk) as can 
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be seen from Eq. ([LI]) . In terms of the canonical partition 
functions, the canonical constraint can be rewritten then 
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If we could relate dZ n /dxjk to Zi(x), . . . , Z„_i(x) then 
the above identity would imply a recursive method for 
computing Z n (x). 

In fact dZ n /dxjk is precisely proportional to a par- 
ticular partition function as can be seen by combining 
Eq. @j and Eq. ((§), yielding 
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So the identity given in Eq. (|12|) is in fact a recursion for 
construction the partition function: 

nZ n (x) = ^ kx jk w kEj Z n - k (x) = ^/cx fc Z„_ fc (x) , 

jk k 

(14) 

with Zq (x) = 1 as the anchor for the recursion. The 
expectation values for njk are simply related to the par- 
tition functions: 



(^jk)n = W k£3 X jk 
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C. Microcanonical partition function 

A similar combinatorial argument can extract the 
microcanonical partition function from the canonical 
or grand canonical partition function. Let us define 
7r(n, E) = {ir jk : kn jk = n, E jfe k ^^ 3 k = E}. Then 



Z n , E (x)= n 



j(n,E) j 
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The partition function can now be developed by recursion 
by noting that 



{^jk)n,E — Xjk 



Z n -k,E-ksj ( x ) 
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Zn—k—m,E — k€j —mei ( x ) /-t r»\ 

= x jk xim ^T~\ , (19) 

and applying the identity V fe k(^jk}n,E = n, yielding 
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with Z ,_e(x) = 5 ^ E - 

A similar situation arises in models of nuclear frag- 
mentation, where there is also two constraints, namely 
the conservation of proton and neutron number, in which 
case [ fl7| it is useful to restrict the ensemble by summing 
over isotopes or isobars. Then the existence of permuta- 
tional Cauchy identities simplified the resulting partition 
functions. In this case the same kind of mechanism also 
applies. We can sum over the energy states and arrive 
at the 7Tfc cluster representation of the states. Or we can 
sum over the cluster sizes and arrive at the rij occupation 
number representation. In each case, the weight retains 
its form and the analysis in the restricted space is iden- 
tical or nearly identical as in the large space. 



III. APPLICATION TO THE HARMONIC 
OSCILLATOR 

Having the partition function now allows us to com- 
pute the thermodynamics in the traditional manner. 
For the case of particles moving in a d dimensional 
isotropic harmonic potential, gj and Sj can be evalu- 
ated directly. The Hamiltonian for the system is given 



by h = ZtiPi/2™ 



T,i=ilf/ 2 where (Qi,Pi) 



is the position and momentum in the ith direction. 
So H = Hi where Hi is the Hamiltonian for a 
one-dimensional oscillator acting on (g;,pi). Of course 
Hitjj n (qi) = (n+l/2)tkij'ip n (qi), and by separation of vari- 
ables Hip nu .„ t n d = (ni H 1- rid + d/2)TkJil> ni ,...,n d - If 

we write this as simply Sj = {j + d/2)huj, we see that 
the jth energy level has a degeneracy gj based on the 
number of ways one can write j as an ordered sum of d 
nonnegative integers, i.e. j = ni + . . . + rid- In combi- 
natorics this is known as a cZ-composition of j, and there 
are gj — (" ,+ ^ _1 ) ways of doing this. Substituting this 

into x k = (1/fc) Y,j 9jW k£j yields 



Xk 
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kU+d/2) = 1 x (21) 



k (1 - x k ) 
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where x = exp{— Tiw/fceT} and we have applied the neg- 
ative binomial identity (1 - x)~ d = J2 n >o ^C^" 1 )- 

A. Canonical thermodynamics 

We can now compute the canonical partition function 
for the harmonic oscillator recursively using Eq. (p^). 
Of interest is the fact that in one dimension the re- 
cursion yields a particularly simple form due to a par- 
tition theorem of Cayley (cf. JTq] , pg. 209), namely: 
Z n , d =i(x) = W 2 /((l - x){l -IP) • • ■ (1 - x n )). This 
partition function was investigated in a different context 
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by A. Z. Mekjian and S. J. Lee p9| . In two and higher 
dimensions, the partition function does not have a simple 
closed form, but can be computed by the above recursion. 

Having obtained the partition function we can compute 
any thermodynamic function of interest by a judicious 
use of partial derivatives. For example, in the canonical 
ensemble, the internal energy U = fcsT 2 (9/9T)y In Z n 
and the specific heat Cy = (dU/dT)v for a group of 
bosons are given by 



U 
C v 



ET dx k 
— — fa), 
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The pressure P = kgT{d / 8V)t^ Z n and the incom- 
pressibility 1/kt — —V(dP/dV)T can be computed in 
a similar fashion. For the harmonic oscillator, the par- 
tial derivatives of Xk with respect to T are readily deter- 
mined. The partial derivatives with respect to V can also 
be determined, once an appropriate volume is specified. 
Since hu> determines the size of the trap, fiui determines 
the volume. Assuming h 2 /2mX 2 w mw 2 A 2 /2 where A is 
a particle's wavelength we see that Tiu oc V~ 2 l d . 

For computational purposes, the above expressions 
are concise and powerful, but for understanding the 
source of the phase transition, one needs to rewrite 
the above expressions in terms of the occupancies of 
the levels. Since Z n = J2n exp{— J2j Tij/kBT} it 
can be shown that U = ^2j( n j)n £ j' ana - that CV — 
Y,jk (( n j n k)n ~ (nj)„(n k ) n ) £j£k/ \k B T 2 ). As the criti- 
cal point is marked by a peak in the specific heat, from 
the above expression its source must be due to large 
fluctuations in the occupation numbers. The expected 
occupation numbers can be inferred from Eq. (|^) and 
Eqs. @ and ©, yielding 



kEj Z n -fc(x) 
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and allow the expected occupation and its fluctuation to 
be computed readily. 



B. Microcanonical thermodynamics 

In the microcanonical ensemble, the temperature is de- 
fined by T = (dE/dS) n y, and this allows us to compute 



the thermodynamic functions, e.g. for the harmonic os- 
cillator 



C v = 



(S n ,E+hu) — S n ,E)\S n ,E ~ S n ,E-hw) 



E-hu 



(26) 



where S n ,E = fog-^n.B- Relating the thermodynamics 
to the occupancies of the levels is more difficult than 
in the canonical case, since the temperature is not an 
independent parameter. 



IV. DISCUSSION AND CONCLUSION 

As an application of the above techniques, we display 
the thermodynamics of particles in a harmonic trap in 
the three ensembles. The specific heat (Fig. 0(a)) and 
the occupation of the ground state (Fig. [|(b)) are sub- 
stantially in agreement in all three ensembles, confirming 
the essential validity of the use of the different ensembles 
even for such small groups of particles. The finite size 
effects are however measurable, and should not be com- 
pletely discounted, especially in smaller systems. 

In summary, we have introduced a combinatorial ap- 
proach to computing thermodynamic functions in sys- 
tems of identical bosons. It is ideally suited to sys- 
tems in which the energy levels are separated by inte- 
gral multiples of some fundamental energy, such as in 
the harmonic oscillator (sj = jtiuj) and the rigid rotator 
(e.j = j(j + l)h 2 /2I). Such techniques can also be applied 
to groups of identical fermions, which has applications to 
the nuclear shell model at finite temperatures. 
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FIG. 1. The specific heat and ground state occupation in three dimensions in the grand canonical, canonical and micro- 
canonical ensembles for a system of 100 particles. The curves become nearly identical for larger groups of particles. 
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